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Abstract. The process of deriving an interatomic potentials represents an attempt to integrate 
out the electronic degrees of freedom from the full quantum description of a condensed matter 
system. In practice it is the derivatives of the interatomic potentials which are used in molecular 
dynamics, as a model for the forces on a system. These forces should be the derivative of 
the free energy of the electronic system, which includes contributions from the entropy of the 
electronic states. This free energy is weakly temperature dependent, and although this can 
be safely neglected in many cases there are some systems where the electronic entropy plays a 
significant role. Here a method is proposed to incorporate electronic entropy in the Sommerfeld 
approximation into empirical potentials. The method is applied as a correction to an existing 
potential for titanium. Thermal properties of the new model are calculated, and a simple 
method for fixing the melting point and solid-solid phase transition temperature for existing 
models fitted to zero temperature data is presented. 



1. Introduction 

Classical molecular dynamics simulations are widely used in many areas of the physical sciences. 
Eliminating the explicit treatment of electronic degrees of freedom brings a computational 
advantage of around six orders of magnitude in the simulation of metallic systems, large enough 
to tackle systems which would be crippled by finite size effects in an electronic structure 
calculation. For applications in microstructure, plasticity and radiation damage where one 
wishes to study the behaviour of the atoms rather than electronic structure, this makes classical 
MD the method of choice. 

The price for this is that the interatomic forces must be represented by the derivative of 
an interatomic potential, which is necessarily a poor approximation to the electronic structure. 
Typically this potential involves some functional form fitted to empirical or ab iratio-calculated 
data. It is useful to distinguish the two distinct sources of error involved in this process: incorrect 
functional form and poor fitting. 

Many MD simulations are still conducted using a pairwise-additive interatomic potential, 
such as Lennard Jones or Coulomb charges. Such pairwise potentials constrain possible values 
of the elastic constants. Most notably, the "Cauchy" relation which relates C12 with Cqq. In a 
pairwise potential, elastic constants are given by the second derivative of the energy with respect 
to strain. Regarding the potential as a function of r 2 rather than r, it is easy to show that for 
any pair potential: 



C 12 = C 66 = ^F"(4)44 



'J 



where rij,Xij,yij are the separation between atoms i and j, and x,y components thereof, where 
i,j run over all atoms and f2 is the volume of the system. Experimentally, this is true only for 
strongly ionic materials. Table 1 shows that elasticity in NaCl can be described by pairwise 
forces, but it is impossible to fit the three independent elastic constants for typical oxides, noble 
gases, water and metals. No amount of fitting can circumvent this mathematical constraint: the 
problem is one of of incorrect functional form. 



Table 1. Elastic constants for selected materials showing the violation of the cauchy relation 
C12 = 6*66 and by implication the impossibility of describing these materials by a pairwise 
additive potential. For water the constants are C13 and C44, the hexagonal equivalent relation 



Material 


C n 


C12 


Cm 


NaCl 


482 


128 


127 


MgO 


291 


96 


152 


Argon 


233 


149 


117 


SiC 


385 


135 


257 


Cu 


168 


121 


75 


Fe 


237 


141 


116 


Water ice 




59 


31 



There are many similar examples where an incorrect functional form renders the process of 
fitting impossible. About 25 years ago there was a strong move away from pairwise forces, 
in particular in metals, in favour of models which describe the local environment such as the 
Embedded atom method[3j, Finnis-Sinclair[2], and Tersoff[4j. Easily justified by elementary 
consideration of the electronic structure, these functional forms allowed enough freedom in 
fitting to avoid the Cauchy constraint and similar relations governing surface and vacancy 
energies. They satisfy the needs of most applications, and even if transferrability is imperfect 
have remained the workhorse of atomistic simulation ever since. 

Some more recent potential developments have simplified the fitting process or made it 
more intuitive. The Modified Embedded Atom Method[5] maps local coordination onto 
spherical harmonics and hence establishes a link to electron orbitals. Two band and magnetic 
potentials [6l El [8] explicitly introduce limited electronic degrees of freedom, and then show how 
to eliminate them analytically, leaving a simple potential of embedded atom form, but justifying 
more complicated parameterisations. 

Ultimately, however, derivation of a functional form for a potential should go back to first 
principles. The overwhelming success of the Kohn-Sham formulation of the density functional 
theory suggests that we should start with the Kohn-Sham Hamiltonian[9j, which in standard 
notation is: 

(1) 

For molecular dynamics we are interested in forces. These can be calculated using the Hellmann- 
Feynman theorem [TO] which tells us that the forces are simply the partial derivative of the Kohn- 
Sham Hamiltonian with respect to ionic positions R{. Inspection of eqUJ shows that only the 
final two, electrostatic, terms contribute to this. 

At finite temperature the electron free energy for a metal changes as the Fermi distribution 
becomes non-singular. The energy increases as states above the Fermi are occupied, while the 



electronic entropy also increases due to partial occupations. According to Sommerfeld theory [TTj 
there is an additional temperature-dependent contribution to the free energy which depends on 
temperature and the density of states at the Fermi energy as 

F som {T) cx T 2 n{E F ) 

In many electronic structure packages this effect is exploited to improve the numerical stability 
of the self consistency loop, by calculating with finite temperature electrons, so-called "Fermi 
Smearing". Typical effective temperatures can run to thousands of Kelvin, and the result is 
adjusted back to zero Kelvin. Nevertheless, the Sommerfeld temperature dependence of the 
electronic free energy is real, and so any potential involving integrating out electronic degrees 
of freedom should itself be temperature dependent. It is worth recalling at this point that the 
Hellmann-Feynman theorem is valid only for an electronic energy which is a variational minimum 
with respect to the parameters describing the electronic structure. For the Fermi-smeared system 
the variational quantity is the electronic free energy (U e i — TS e i), so the Hellman-Feynman 
forces are the derivative of this. The physical assumption underlying this is similar, though not 
identical, to the Born-Oppenheimer approximation: the ionic motion should be slow enough for 
the electrons to relax into their equilibrium distribution. 

The temperature dependence of the Sommerfeld electronic free energy for Ti is shown in FigQ] 
alongside the associated density of states. These figures were calculated using the CASTEP[T3] 
program with settings as given in previous work [I], using Fermi-Dirac smearing with the 
smearing width corresponding to the quoted temperature. In all cases the atoms were located 
on their ideal lattice sites. The density of states was calculated using the castepdos program, 
which interpolates and integrates the density of states using the octahedron method between 
the explicitly calculated k-points.[14j At 0.05eV (about 600K) the Sommerfeld contribution in 
bcc is 34m «V nnmrmrprl with 7meV fnr hen This is about 20% of the free energy difference at 
OK. 
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2. Sommerfeld Potentials 

In order to include Sommerfeld effects in molecular dynamics, the potential must be written in 
the form of a sum of energies per atom j in an arrangement of atoms i in positions {r;} : 

C/({ ri }) = tfbfa, {>i}) + F sam {v h { ri }) (2) 

j 



The vast majority of interatomic potentials are fitted to (extrapolated) zero-temperature 
experimental data and/or ab initio data calculated on the Born-Oppenheimer surface. Thus 
in many cases good parameterisations for Uq already exist. Here we consider the second term in 
equation [21 which is proportional to T 2 . The proportionality constant can be treated as a fitting 
parameter (At). The challenge is to obtain a sensible measure for the local density of states at 
the Fermi energy rij(Ef) in terms of the atomic positions without performing a full electronic 
structure calculation. The most convenient approach is to write it as a pairwise function: 

^(£/) = E/(i r i- r ji) 

i 

We could stop at this point, and simply take / as an arbitrary function to be fitted, however it 
is worth discussing how it might look. 

2.1. Heuristics for f(r) in transition metals 

f(r) attempts to measure the density of states at the Fermi level projected onto an individual 
atom. This quantity can be readily calculated for various crystal structures (e.g. Fig 2), and 
may vary quite widely between them. If we consider a canonical d-band model, the density of 
states is very different between close-packed and bcc structures. The DoS can be defined in 
terms of moments, in particular the third and fourth moment determine the skew and kurtosis. 
The moments theorem [16] relates these to a real-space picture involving closed paths of near 
neighbour hops. This defines the key difference between close packing (many loops of three 
hops) and body centred cubic (no three-hop loops- many four loops). Open structures seldom 
feature in transition metals and their alloys, so we can neglect these. 

Although counting paths of hops cannot be done directly by a pair potential, it suggests a 
related local measure. For practical purposes the skew/kurtosis relation can be approximated 
by using a function which is either sharply peaked around first neighbour (favours three-loop 
equilateral triangles) or a broader function which favours four-loops. Which is energetically 
favoured for a given material depends on the band filling. 

In most existing EAM and Finnis Sinclair potentials this narrow vs broad minimum in the 
potential determines the stability of close packed vs bcc materials. In the moments picture, 
the longer range acts as a proxy for measuring four-membered rings. Use of exponential fitting 
functions, as in early EAMs, tends to give a narrow minimum; polynomial functions, as used in 
Finnis Sinclair potentials, give a broader minimum. This, rather than any deep physics, explains 
their early application to fee and bcc respectively. 

2.2. Application to Titanium 

For titanium the Fermi level falls at a minimum of the hep density of states, and a maximium 
of the bcc density of states. Thus we can expect that electronic entropy contributions will be 
different and significant in determining the phase transition temperature. Indeed, electronic 
entropy was shown to provide almost half the excess entropy of bcc Zr[12], a material very 
similar to Ti. The density of states is for valence electrons, so f(r) should have significant value 
only outside the repulsive core. It is also desirable that its derivatives be continuous and the 
number of adjustable parameters is minimised. Consequently we write the function as: 

f(r) = X 2 (l-X) 2 0<X<1; X = (r-r )/d 

Fsom = X] A T T 2 X 2 (1 - X) 2 

ij 

this gives a three parameter model, defining the range of the interaction from ro to r Q + d and 
the strength At 



Most published potentials are fitted to zero temperature data, so this Sommerfeld term 
can simply be added. Moreover, the parameterisation is relatively straightforward: we have 
calculated the temperature-dependent contribution to the free energy for fee and hep Ti. These 
two pieces of data are sufficient to parameterise the model. 

For Ti we consider a potential introduced in 1992[21j which is known to have too low a basal 
stacking fault but otherwise reproduces the behaviour of hep titanium reasonably well. From 
Fig. UJ above we see that at low energies (up to 1000K) the Sommerfeld contribution to bec is 
about 4.7 times that of hep. This suggests an approximate solution r Q = 2.8AA, d = 1.46A, 
A T = -7.5 x 10- 7 meV/K 2 . 

2.3. Which temperature to use? 

In order to use the temperature-dependent potential, it is necessary to specify the electronic 
temperature. In practice, most MD simulations are relatively small and one can assume that 
the electrons are able to reach thermal equilibrium throughout the region. 

Thus for an NVT/NPT calculation run using a thermostat, one can simply use the thermostat 
temperature to fix a constant potential throughout the run. For NVE one can use the 
instantaneous global temperature of the simulation for all atoms, so although the potential 
may vary in time, it is constant across space. 

3. Tests of the new potentials 

To test the thermal effects we have carried out melting point calculations from bec and hep Ti, 
and thermal expansion calculations. 

3.1. Melting Point Calculations 

We determined melting point using the coexistence method [18], using the MOLDY program in 
the NPE ensemble [1 7] . By using intermediate sized systems, 16000 and 13824 atoms for hep and 
bec respectively, it is possible to suppress the bec-hep transition and calculate melting points 
from both hep and bec phases. The original potential[21J gives a melting point of 1395 ± 10K 
(hep) and 1790 ± 10 K (bec), the difference indicating the stability of bec at high temperature. 
This contrasts with the OK energies of -4.853eV (hep) and -4.807eV (bec) which indicate low 
temperature stability of hep. The transition at intermediate temperature is via a soft phonon 
analogous to zirconium [121 ITS'] , and care must be taken to ensure that the solid phase does 
not transform during the simulation, which we did by monitoring local coordination using the 
BallViewer code 19 . 

Experimentally, the melting temperature of Ti is 1941K, and the bec-hep transition is at 
1155 K. The original potential therefore gives too low a melting temperature. The effect of the 
Sommerfeld correction on the melting is small (1315K and 1800K), slightly stabilising bec and 
destabilising hep. 

Metals melt at high temperature because the liquid phase has higher entropy than the solid. 
Thus the more of the phase space the system can reach without too high potential energy, the 
lower will be the melting point. Since a melt samples interatomic distances shorter than typically 
realised in the solid, softening the potential inside nearest neighbour distances can be expected 
to favour the liquid. In the present potential, the ag spline parameter controls the short range 
repulsion. The 1992 potential was not fitted to any data for short-ranged interactions, and in 
Fig. [3] we show the evolution of the melting temperature as a function of the parameter a§: the 
short ranged spline starting close to nearest neighbour separations (see Appendix A) . The effect 
on the melting point is very pronounced, while the OK properties are essentially unaffected. 

The original potential had a$ = 0.494eV/A 3 : the experimental melting point (1941K) and ab 
initio hep-bec difference (90meV) were not fitted in 1992, and to do so requires a significantly 
stiffer short-ranged potential. 
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Figure 3. Effect of short ranged repulsion of melting point and zero-temperature hep-bec 
stability. See Appendix A for definition of a$ 

3.2. Thermal Expansion 

The thermal expansion is not typically fitted in empirical potentials, and in consequence is often 
significantly wrong. The Sommerfeld correction introduces additional temperature-dependent 
anharmonicity into the potential, which affects the thermal expansion. For this potential the 
thermal expansion is close to the experimental data but highly temperature dependent, especially 
at high temperature where the softening around second neighbours leads to high values. 

3.3. Phase Transformation Calculations 

Whatever value is used, the Sommerfeld correction will be insufficient to reduce the enthalpy 
of bec below hep: the bec structure is always stabilised by phonon entropy. The constraints 
imposed by the crystal symmetry mean that it is not possible to calculate the phase transition 
temperature between hep and bec by coexistence. Rather, it requires a full free energy 
calculation [2(31 121] which can be performed by integrating using the Gibbs-Helmholtz equation: 



The right hand side of this equation can be calculated from a single-phase NPT molecular 
dynamics calculation (see e.g. Appendix B), varying the temperature slowly by incrementing 
the required temperature of a Nose thermostat. In practice we use this equation with the two 
melting points and assume that H/T 2 varies linearly with T. The stabilising effect on bec lowers 
the transition temperature considerably, the actual value depending on the chosen value of a§. 

Figure 0] shows the effect of the Sommerfeld potential in a dynamics calculation of the phase 
transition. Although this is only a lower bound for the thermodynamic transition temperature, 
it does illustrate the increased stability range of bec. 

4. Conclusions 

The thermal excitation of electronic degrees of freedom introduces a temperature dependence 
into the electronic free energy. If the electronic degrees of freedom are integrated out to make an 
interatomic potential, this temperature dependence should remain. Where the excitation arises 
simply from the Fermi-Dirac distribution via Sommerfeld theory, this temperature dependence 
scales as T 2 . The magnitude of the effect has been calculated for Ti using density functional 
theory, and found to give a contribution of a few tens of meV to the bec-hep free energy difference. 
This in turn is enough to lower the calculated hep-bec phase transition temperature significantly, 
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Figure 4. Plot of cohesive energy per atom vs temperature for MD simulation on heating 
with and without Sommerfeld correction. Main figure shows linear increase in energy with 
temperature, as expected from virial theorem (original potential, black line). Corrected potential 
(red) shows deviation from virial behaviour as electronic degrees of freedom are excited at high 
T. Both cases have a discontinuity as the cooling bcc transforms to twinned hep, the correction 
stabilises bcc, dropping the observed transformation from 1000K to 650K. Inset shows pair 
potential V(r) at OK and 1250K . 



but has less effect on the melting temperature. Neither effect is large compared to the errors of 
typical EAM-type potentials which do not include melting temperature as a fit parameter. 

By contrast, the melting temperature has strong correlation with the short ranged part of 
the potential, inside the normal separations found and fitted at OK. By adjusting only the short- 
range terms it is possible to refit the melting point of a potential fitted to elastic moduli, cohesive 
energy, vacancy and surface formation etc. without disturbing those properties. This suggests 
a simple way to improve the high temperature phase performance of existing potentials without 
major redevelopment [26| 127} 128] . Note that energies of self interstitials[2j [3l [7] and interstitial 
impurities [14^ [29] may be affected. 

In addition to Sommerfeld correction, a strong temperature dependence of the electron free 
energy may be found in magnetic materials such as iron|30[ [31] [32], where the fee phase 
is stabilised by paramagnetic free energy. This effect could be captured by a temperature- 
dependent potential, although the excitation of all magnetic modes means the temperature 
dependence would be more complex, requiring an additional many-body term rather than a 
pairwise one. This will be the subject of future work. 

5. Appendix A 

The potential has a functional form for the energy of the i th atom 



U i = Y,H H ^k- R i 3 ) a k(rk-Ri 3 f (3) 



3 k 



Y,A T T 2 X 2 (1 - X) 2 H(r + d- Ri^HiRij - r Q ) (4) 



3 



- T,T, H ( R k- R ij)MRk-R l3 ) 3 (5) 

V 3 k 

with X as defined above, Rij the separation between i th and j th nearest atoms, H the Heaviside 
step function and the other parameters given inTable Al below. 
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ak 


-0.785715 


1.110966 


-0.299450 


-0.143061 


1.025368 


0.494293 


rk 


5.09113 


5.00767 


4.673828 


3.964408 


3.338449 


2.9508 


A k 


0.547614 


-0.551266 










Rk 


5.09113 


4.381714 












A T = -7.5x10"' 


r =2.84 


(2=1.4 









Table 2. Table Al. Parameters for the Ti potential, from ref.[2T], converted into units of eV 
and Angstroms, and for the temperature-dependent correction. Note that a§ can be varied to 
fit th 
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Figure 5. Plot of enthalpy vs temperature showing failure of thermodynamic integration 
due to phase transition from bcc to deformed hep on cooling. Upper line (black) is bec at 
high temperature, transforming to twinned hep on cooling. Reheating gave hysteresis in the 
retransformation. Note the different transuition point compared with fig 4. Lower line (red) 
is hep on heating from OK: recooling produced the same curve. NPT simulations were run 
with a cooling rate of lK/ps and a lfs timestep, with a Nose relaxation time parameter of 
lOOfs. Phases were initialised in a supercell with periodic boundaries compatible with the 
phase in question: although the transformation could happen in principle in Par rinello- Rahman 
dynamics, in practice the transformation is slow enough to gather phonon statistics. 

6. Appendix B 

One potential drawback of the thermodynamic integration is that the system will undergo a 
phase transformation as the temperature is varied, rather than remaining in the metastablc 
phase. When this occurs, the latent heat causes a discontinuity in a graph of enthalpy vs 
temperature. Figj5j shows an example of this in cooling of the bcc phase. Some authors 
have maintained that if the transformation occurs quickly enough, the temperature of the 
discontinuity would be the transformation temperature. However, better statistics than are 
available here are needed to ensure that this is so [25]. Moreover, Fig [5j shows that the 
transformation is not to the pure hep phase, and BallViewer analysis |19j shows that the low 
temperature phase is a twinned hep. 
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